Improved Local Lattice Approach for Coulombic Simulations 
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An improved approach to the simulation of strongly fluctuating Coulomb gases, based on a local 
lattice technique introduced by Maggs and Rossetto 1], is described and then tested in a problem 
of biophysical interest. The low acceptance rates for charged particle moves in regimes of physical 
interest are increased to a serviceable level by use of a coupled particle-field update procedure in 
the new method. Sensitivity of the results to lattice discretization effects is also studied using 
asymmetric lattices. 



I. INTRODUCTION 

The simulation of systems containing a large number 
of mobile charged entities, in which long range electro- 
static forces play a central dynamical role, is of critical 
importance in modern chemical physics and biophysics. 
In many cases, the computational load is dominated by 
the evaluation of the electrostatic energy of the system, 
where the long-range character of the Coulomb interac- 
tion greatly complicates the development of efhcient algo- 
rithms that scale with system size in a way that permits 
study of systems of biophysical interest. Developments in 
supercomputing technology aimed at large-scale biophys- 
ical simulations, such as the IBM BlueGene project [31 , 
where massively parallel assemblies of processor nodes 
are coupled via a three-dimensional toroidal topology, 
suggest that algorithms based on a local energy func- 
tional will be much more efficiently executed on the next 
generation of high-end computing platforms than those 
involving long-range nonlocal effects. 

Recently, Maggs and collaborators Qj have suggested 
an ingenious procedure for removing the nonlocal (long- 
ranged) Coulomb term in equilibrium simulations of 
Coulomb gases. By using a completely local Hamilto- 
nian for a system of mobile charged particles interact- 
ing with the electrostatic field, one avoids the unpleasant 
scaling characteristics of conventional Coulomb gas sim- 
ulations. Unfortunately (as pointed out by these authors 
themselves Q), the algorithm they propose runs into se- 
rious acceptance problems in regions of physical interest 
(basically, for strongly fluctuating systems). In this pa- 
per, we study the origin of these acceptance problems and 
propose an improved algorithm that allows useful simula- 
tions of strongly fluctuating systems in which mean-field 
(or Poisson-Boltzmann) methods break down. 

In Section 2 we briefly review the original technique 
of Maggs et al, and explain the origin of the acceptance 
difficulty for charged particle moves. In Section 3 we ex- 
plain the modified update procedure designed to cure, 
or at least ameliorate, the acceptance problem for par- 
ticle moves. In brief, the crucial point is to implement 
a coupled particle-field update in which the electrostatic 
field is allowed to readjust itself in tandem with charged 



particle moves in response to the changed electrostatic 
environment. In Section 4, the improved procedure is 
tested in detail on a system that has been extensively 
studied in the literature 0, IE IE ■ the osmotic pres- 
sure of charged plates (or membranes) separated by an 
ionic fluid. Finally, in Section 5 we briefly summarize our 
conclusions. 



II. LOCAL LATTICE HAMILTONIANS FOR 
COULOMB GAS PROBLEMS 

The difficulties incurred by the nonlocal nature of 
the Coulomb interaction in realistic simulations of large 
systems (for example, for large biomolecular systems) 
are well known: the computational cost increases as 
the square of the number of charged constituents, and 
although various techniques (Ewald summation, fast 
Fourier transforms etc. [Sj) can be employed to improve 
this scaling, the resulting complications in the algorithm 
often mean that the computation of the electrostatic en- 
ergy still consumes essentially all of the computational 
effort, greatly limiting the size of the systems and (in the 
case of molecular dynamics simulations) the time frames 
over which the simulations can be extended. These tech- 
niques also have difficulties modeling a non-uniform di- 
electric constant, which is an important feature of many 
biophysical systems|E 0| as the dielectric constant 
in proteins is '--^ 2 — 8 while the dielectric constant of wa- 
ter is ~ 80. In the case of systems at equilibrium, it has 
been known for some time y,Uj,Lj,Uj that the nonlo- 
cal Coulomb interaction can be replaced by a completely 
local interaction via a Hubbard-Stratonovich transforma- 
tion, yielding a path integral formalism that connects 
naturally with the Poisson-Boltzmann mean-field theory. 
Unfortunately, for strongly fluctuating systems pertur- 
bation theory (saddle-point expansions) breaks down in 
this approach, and a direct numerical simulation is ob- 
structed by a severe sign problem. 

Recently, Maggs and collaborators ij have proposed 
an alternative, purely local approach to the simulation 
of charged condensed systems. They exploit the fact 
that the nonlocality of the Coulomb interaction is a con- 
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sequence of a particular choice of gauge for describing 
the electromagnetic field, whereas the physically relevant 
quantity -the electrostatic energy of the system- must 
clearly be a gauge-invariant object. They propose that 
the electromagnetic field be simulated in terms of gauge- 
invariant objects (specifically, the electric field), repre- 
sented on a discrete spatial lattice. In this respect, the 
method proposed is essentially the same as that employed 
for over 20 years by elementary particle theorists attack- 
ing the problem of strong interactions with the technique 
of lattice quantum chromodynamics. For a review see 
Ref. 13 The main distinction here is that the gauge 
theory involved is the abelian one of Maxwellian elec- 
trodynamics, magnetic effects are not relevant, and the 
formulation used is a noncompact one (i.e. the electric 
field variables take unbounded values). 

Let us briefly recall the salient points of the formalism 
of Maggs et al . The canonical partition function for 
a set of mobile charges at locations Vi at inverse tem- 
perature /3 in a medium of dielectric constant e may be 
written 

Z= !\{ dn VE{r) \{5(^-E~ -p(f)) e-fe / 
•' i=i f ^ ^ ' 

(1) 

where the charge density p(r) is shorthand for 

p(f) = ^e,5(f-rl) (2) 

i 

The delta function constraint in Eq. 2] enforces Gauss' 
Law, so that the electric fields integrated over correspond 
to the particle locations specified through the density 
function p. The formulation is manifestly local, as both 
the energy functional and Gauss' Law constraint are so. 
There is no requirement that the electric fields integrated 
over be irrotational, and in fact they are not; as shown 
by Maggs et al.0, the transverse part of the electric field 
simply decouples from the particle sector and contributes 
an irrelevant overall prefactor to Z. Because this for- 
malism is local, it is easily extended to model physical 
systems with a non- uniform dielectric constant 16]. 

The functional integral over electric field in Eq. ^ can 
be given a precise definition by introducing a spatial cu- 
bical lattice, which we shall for the time being take to be 
a grid of 1? points, with lattice spacing a (in all direc- 
tions: the modifications needed in case of an asymmetric 
lattice are discussed below) and periodic wrap-around 
boundary conditions in all three spatial directions. The 
charges e^, i = 1, ..iV are assumed to be integer multiples 
of a basic unit of charge, = z^e, zi integer, and re- 
side on the sites of the lattice. The component Ef^{fi) of 
electric field in direction /i at lattice site n is associated 
with a real-valued field Ei on the oriented link I from n 
to fi + ji. Discretizing the 3-dimensional integral for the 
electrostatic energy in the obvious way, we find 



The implementation of the simulation is simplified by 
introducing dimensionless variables to the greatest extent 
possible, so we define Ei = ^^Ei and a rescaled inverse 

temperature (3 = in terms of which the energy 

becomes 

H-^j:^f (4) 

while the Gauss' Law constraint takes the simple form 

E^' = ^^ (5) 
I 

for the sum of outgoing link fields from any site contain- 
ing a charged particle of charge z^e. 

The simulation of the system defined by the energy 
function in Eq. ^and the constraint in Eq.jSjcan in prin- 
ciple be accomplished by the following algorithm: 

1. Pick starting lattice locations (possibly randomly) 
for the N particles of charge Zi,i = 1,..N. Then 
solve Gauss' Law for these fixed charge locations 
to obtain a starting configuration of electric link 
field variables satisfying the Gauss constraint. This 
can easily be done by standard numerical relaxation 
methods l^. 

2. Update the electric fields by shifting all link vari- 
ables along a complete set of independent closed 
paths by constant shifts, using either Metropolis 
or heat bath procedures to accept /reject proposed 
shifts. The simplest version of this is simply to 
consider all plaquettes (unit squares) on the lattice, 
shifting the 4 link fields ordered around the plaque- 
tte by the same random amount a, the range of a 
set so that there is a reasonable acceptance rate for 
the move. Such a shift clearly maintains the Gauss' 
Law constraint. 

3. Update particle locations by visiting in turn every 
site n containing a charged particle of charge Zi. 
A particle move to the neighboring site n + jl in a 
random direction fi is then considered, where the 
particle move is accompanied with a shift of the 
electric field Ei on the link I = {n ^ n + fi) 

Ei^Ei- z, (6) 

in order to maintain the constraint in Eq. Here 
also one can employ either Metropolis or heat-bath 
accept/reject procedures. 

The inclusion of additional force fields, for example soft 
or hard exclusion potentials modeling a finite size for the 
particles, is, in principle, completely straightforward in 
this framework. When particles are packed closely to- 
gether, or the potential changes rapidly over the scale 
of a lattice spacing, then it is important to verify that 



the observed phenomena are not distorted by lattice dis- 
cretization effects. It is useful to be able to study the 
effects of lattice discretization in such situations by in- 
troducing asymmetric lattices, in which the lattice spac- 
ing in the various directions differs. As a specific exam- 
ple, consider a situation in which we may desire a finer 
discretization in the x-direction, relative to the y and z 



directions, ax < ay 



One readily verifies that 



with the choice of the dimensionless variables 
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(7) 



where La is the set of links in the a direction and, as 
before. 
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the energy function becomes 




^Ef + y: -Ef + E -^f 



(8) 



(9) 



while the Gauss' Law constraint retains its original form 
given in Eq. 

Unfortunately, despite the appealing simplicity of the 
simulation procedure outlined above, in physically real- 
istic situations involving strongly charged systems the 
method proves impractical, for reasons we now ex- 
plain. The dimensionless inverse temperature variable 
/3 is typically large compared to unity (in the charged 
plate/membrane problem considered in Section 4, the 
value is 87.1), so that typical values for the electric field 
link variables are small compared to unity. On the other 
hand, executing a particle move across a link via Eq. |^ 
shifts the electric field variable on that link by an integer, 
and this generally leads to an unacceptable energy cost 
(on the order of (3). In the univalent case {zi = ±1), ac- 
ceptance rates for particle moves are of the order of 10"'*, 
while for divalent ions {zi = ±2) the acceptance rate is 
at best of order 10~^. Thus the unmodified procedure 
of Maggs et al. is clearly not a practical approach in sit- 
uations approximating real biophysical systems. In the 
next section, we discuss a modified simulation algorithm 
in which this problem is ameliorated to an acceptable 
level. 



III. SOLVING THE PARTICLE MOVE 
PROBLEM: A COUPLED UPDATE PROCEDURE 

The problem of very inefficient particle moves men- 
tioned in the preceding section needs to be resolved 
before the local Hamiltonian method can be applied 
fruitfully to realistic problems with strongly fluctuating 
Coulomb gases. Recall that the Hamiltonian, as a func- 
tion of the electric field variables Ej defined on the links 




FIG. 1: Local field environment for coupled particle move 
updates 



I of the lattice takes the form 



(10) 



where the dimensionless inverse-temperature variable, /3, 
is quite large for the systems that we are interested in 
studying (in the range of 50-100). As discussed above, the 
vast majority of particle moves with such a Hamiltonian 
have a high energy cost leading to an unacceptably low 
acceptance rate. In this section we will show that this 
problem can be substantially ameliorated- though not 
completely eliminated- by a coupled update procedure in 
which electric field values on all the plaquettes containing 
the link along which the particle move is attempted are 
simultaneously adjusted to reflect the changed electrical 
environment resulting from the particle move. 

As an example of a simple procedure that can consid- 
erably improve the acceptance rate for particle moves, 
consider the situation illustrated in Fig. ^ Here we are 
considering the move of a unit charge particle from the 
beginning (bottom) site to the end (top) site of the cen- 
tral link associated with field variable Eq. In conjunction 
with the particle move, we consider simultaneous electric 
field updates corresponding to plaquette variable shifts 
ai,a2,a3,a4 on the four plaquettes containing the link 
Eq, as indicated in the figure. Such a combined move 
changes the energy associated with the illustrated region 
from 
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In practice one finds that the electric field variables equi- 
librate to values which are small compared to unity: in 
the approximation where we simply set .E; = in Equa- 
tions and El the energy cost of the combined move 
becomes 

= § - + + «2 + «i + «4) j 

(13) 

Minimizing Eg. 1131 with respect to the a^, we find that 
the choice = 1/7 gives the minimum energy cost 



(14) 



as opposed to the cost /3/2 if the particle move is un- 
accompanied by any readjustment of nearby link fields. 
As we shall see in the next section, if f3 is large, this is 
enough to increase the acceptance rate for particle moves 
to a level where configurations can be decorrelated at an 
acceptable rate. Thus, a quick and easily implementable 
improvement of the basic algorithm can be obtained by 
a Metropolis accept/reject step in which the choices for 
a particle move on a chosen link are (a) do nothing (to 
particle or fields), or (b) perform the combined update 
in which the particle is transferred to the end site of the 
link and the fields around the four intersecting plaquettes 
are shifted by 1/7. It is clear that the energy cost can 
be further reduced by allowing readjustments of plaque- 
ttes adjacent to those depicted in Fig.^ In particular, if 
the link Eq in Fig. corresponds to the z direction, then 
adding the remaining xz and yz plaquettes that contain 
the links E2, E^, E^ and En, and performing the rele- 
vant minimization, one finds however that the reduction 
in energy cost is only about 10%, with a considerable 
complication in the algorithm. In this paper we have 
chosen to implement only the simplest (most local) ver- 
sion of a coupled move-field update, corresponding to the 
situation in Fig. 

A more general procedure, in which a heat-bath update 
on the combined (particle move) -I- (field update) space 
provides a complete local decorrelation between adjacent 
Monte Carlo configurations can easily be derived as fol- 
lows. We remind the reader that, in a heat-bath Monte 
Carlo update, parameters are introduced to characterize 
a subspace of the configuration space in the neighborhood 
of the starting configuration, the dependence of the full 
Boltzmann weight of the theory on these parameters is 
extracted (from the full Hamiltonian) , and then new val- 
ues for these parameters are then chosen (independent of 
the original configuration) on the basis of this Boltzmann 
weight. In the situation considered here, the parameter 



space consists of a single discrete particle move variable 
m = 0, 1 (with m = corresponding to no move, m — 1 
to a move along a specified link), and four continuous 
plaquette shift variables ai,i = 1, ..4. For the indicated 
environment of the central link Eq in Fig. ^ the rele- 
vant part of the Hamiltonian, as a function of m and the 
continuous plaquette update variables a^, becomes 



7i(m, ai) 



■2|(^o -"iz)^ 

3 6 
Y,{Ei + aif +Y.^Ei + a2f 



(15) 
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1=4 
12 

+ Y.^Ei+ a^f + Y.^Ei+ aif 

1=7 l=W ) 

where we have introduced a variable (integer) valence z 
to take care of the case (needed in the simulations of 
Section 4) of multivalent ions. In order to implement 
a heat-bath procedure for this energy function, we need 
to generate values for the quintet (m, ai, 02, a^, 04) dis- 
tributed according to the Boltzmann weight e~^(™'"'). 
Fortunately, a complete analytic solution to this prob- 
lem can easily be derived. First, we note that the energy 
function in Eq. 1151 can be reexpressed 



H{m, ai) = {Eq - mz^ + ^ aiMijUj (16) 

12 



+2Y,Ka, + Y,E^i 



where 



1=1 



Xi — Pi — mz 

Pi = EQ + E1+E2 + E3 

P2 = E0 + E4 + E5 + E6 
etc 



(17) 
(18) 
(19) 
(20) 



i.e. the variables Pi, i = 1,2,3,4 are just the plaque- 
tte fields obtained by summing the electric link variables 
around each of the four plaquettes containing the central 
link of Fig. ^ and the 4x4 matrix My 



is easily found to have inverse 



= 1(5, - 1) 



(21) 



(22) 



Completing the square in Eq. 1161 we find that H takes 
the form 



I ij 



(23) 
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ij 1=1 ) 

^ a,+Y,Mr:^X, (24) 
j 

The dependence of the local energy on the discrete move 
variable m arises from the first and third terms in Eq. 1231 
and the corresponding Boltzmann weight determining 
the relative probability of a particle move (m = 1) versus 
no move (m = 0) is therefore 




where we have used the fact that the move variable m — 
0, 1 so that — m. A heat bath update of the variable 
m is therefore trivial to implement. 

The continuous ai variables can be generated easily 
from the Gaussian distribution of the 6- The eigenval- 
ues of Mij are easily found (they are 7,3,3,3), as are the 
eigenvectors, whence we find that the contribution of the 
second term in Eq. 1231 to the Boltzmann weight can be 
rewritten 

exp 1 {77^1 + 377I + 3%2 + 377I) ^ (26) 

where 



Vi = 


^(6 + 


6 + 6 


+ 6) 


(27) 


m = 




-6) 




(28) 


V3 = 




-6) 




(29) 


774 = 


^(6 + 


6-6 


-6) 


(30) 



The heat-bath procedure for the plaquette shifts a; there- 
fore amounts to generating the independent Gaussian dis- 
tributed variables r/i according to the weight (|26|l . whence 
the ai can be reconstructed via Equations 1271301 and 

= 6-^A, + l^A, (32) 
i 

To summarize, the algorithm for a coupled parti- 
cle/field heat-bath update is implemented as follows: 

1. Calculate the plaquette sums Pi,i = 1,2,3,4 
fEauations I18l20|l for the four plaquettes interfac- 
ing the link along which we desire to move the par- 
ticle. 

2. Choose the move variable m = 0, 1 with weight 
given by Eq. 1251 



3. Generate independent Gaussian variables r]i,i = 
1, 2, 3, 4 according to Eq. ^ 

4. Solve Equations EEni for the ^„i^l,2, 3, 4. 

5. Compute from Ea. 1171 and use Eq. |221to obtain 
the desired plaquette shifts ai,i — 1,2,3,4, which 
are then used to update the electric fields Ei,l — 
0, ..12 as indicated in Eq. ^l(see Fig.nj. 

6. If TO = 1 then move the particle across the con- 
sidered link while updating the electric field on the 
link according to Eq. |H1 

Recently the problem of low acceptance rates for par- 
ticle moves was noted by Maggs et al. in Ref. 0. They 
present an alternative solution to the problem where each 
charge, instead of residing on a single lattice site, is bro- 
ken into pieces and resides on the lattice sites in an nxnxn 
cube. In order to move a particle, all of the pieces of the 
particle must be moved in unison. They have shown that 
the inverse temperature that they are able to simulate ef- 
ficiently grows like using this method. The advantage 
of this method is that it is effective at increasing the ac- 
ceptance rate and that the size of the cube can be chosen 
to give the desired acceptance rate. The disadvantage of 
this method is that the charges are spread out so that, 
for systems that are sensitive to the spatial location of 
the charges, the lattice must be made finer by a factor of 
n in every direction to obtain the same charge locality as 
the lattice with unbroken particles. Using the methods 
discussed in this work, the charges remain on a single 
lattice site so there are no difficulties arising from the 
breakup of the ions onto different lattice sites. 

IV. APPLICATIONS: STRONGLY 
FLUCTUATING FIELDS BETWEEN CHARGED 
PLATES /MEMBRANES 

To test these algorithms on a strongly charged sys- 
tem where correlation effects play a major role, we have 
considered a system of charged conducting plates with 
ions between the plates. While the system is electri- 
cally neutral, there is an osmotic pressure between the 
plates that depends on the electrostatic interaction be- 
tween the particles and on the correlations between the 
particles. This system has been extensively studied both 
theoretically^, |^ and numerically^, Q and is known 
to be a strongly fluctuating system within the param- 
eter ranges in which we are interested. A convenient 
criterionQ for a strongly fluctuating system is that the 
Bjerrum length, £b, times the square of the ion valence, 
z, be smaller then the Gouy-Chapman length, /i. In the 
systems we are considering z^is/fJ' is as large as 33. We 
have considered both divalent and univalent ions, as in 
previous work it was seen that there is a repulsive pres- 
sure in the univalent case and an attractive pressure in 
the divalent case0, 0, 0- -"-^ shown in Ref. Q| that the 
Poisson-Boltzmann calculations of the osmotic pressure 
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in the divalent case break down and cannot even predict 
the sign of the osmotic pressure. 

Our basic system consists of a 50x50x50 lattice with 
periodic boundary conditions in all three dimensions. 
Positive charges are free to move on two fixed plates sep- 
arated in the x direction which extend the entire extent 
of the lattice in the y and z directions, while the region 
between the plates contains mobile counterions ensur- 
ing overall neutrality. The periodicity in the x direc- 
tion is not critical, as quantities observed are insensitive 
to field fluctuations far outside the plates. We choose 
a lattice spacing of 1 A, so that we can study plate 
separations in the range of interest. Using the dielec- 
tric constant of water (e — 80.0) and room temperature 
(T = iQQK) gives a dimensionless inverse temperature 
/3 = 87.1, too large to effectively simulate with simple 
particle moves that do not adjust the electric field on 
neighboring plaquettes. We placed 34 positive univalent 
charges on each of the plates to give a surface-charge den- 
sity of 0.2176 Cm~^, approximately that used in Refer- 
ences 0, 0. These charges are allowed to move during 
the simulation, but are not allowed to leave the plate. 
To make the system electrically neutral, 68 negatively 
charged ions are placed between the plates in the univa- 
lent case, and 34 negatively charged ions in the divalent 
case. Two ions are forbidden from being on the same 
lattice site. The charges on the plates are initially ran- 
domly distributed on the plates, and the ions between the 
plates are initially distributed with half of the ions on the 
closest allowed plane to the right plate, and half on the 
closest allowed plane to the left plate. All runs are com- 
posed of 5,000 Monte Carlo equilibration steps followed 
by 20,000 measurement steps. Each Monte Carlo step is 
composed of a coupled Metropolis update of the electric 
field around each plaquette, (200 x Number of charges 
on the plates) attempted moves of a particle on the plate 
chosen at random, and (20000 x Number of charges in 
solution) attempted moves of a particle in solution. As 
pointed out by Maggs et al.[3|j ^ global update of the 
electric field is also included to ensure rigorous ergodicity. 

To investigate the importance of the mobility of the 
charges on the plates we have also performed a set of 
simulations with the positive charges on plates fixed at 
a random initial distribution. There were no qualitative 
differences between the results of these simulations and 
the results of the simulations with mobile ions on the 
plates that we present here. 

To investigate the errors due to lattice effects, we have 
also studied an asymmetric lattice where the lattice spac- 
ing is a factor of two smaller in the dimension separat- 
ing the plates. This is a 100x50x50 asymmetric lattice 
with a lattice spacing of 0.5 A in the x direction so that 
the total volume of the system remains constant. Two 
ions are again forbidden from being on the same lattice 
site. In the asymmetric case this corresponds to a differ- 
ent hard sphere interaction between the ions than in the 
symmetric lattice, but these differences are in practice 
unimportant, as the ions are so sparsely distributed that 



collision between ions are rare. 

The ions between the plates will naturally accumulate 
on the planes of lattice sites close to the plates. As the 
electric potential changes rapidly in this region, the re- 
sults of our simulation will depend on the details of the 
discretization in this region. As the discreteness of the 
lattice has the largest effect in the region close to the 
plates, we have chosen to forbid the ions from coming 
within 1 A of the plates. This will soften slightly the po- 
tential seen by the ions. On the symmetric lattice, we do 
not allow ions on the planes of lattice sites closest to the 
plates. On the asymmetric lattice, ions are not allowed 
on the 2 planes of lattice sites closest to the plates. 

We are primarily interested in observing the osmotic 
pressure between the plates as we change the separation, 
both for the univalent and divalent ions in solution. As 
derived in Ref. 0, the osmotic pressure can be calculated 
using the expression 

Posm = kT C(0) + F^^/area, (33) 

where C (0) is the ion concentration at the mid-plane and 
i^/^ is the average electrostatic force between the left 
half of the system and the right half of the system. In 
the continuum, this force could be written as 

^ A B 

^^"^ qnqrnAx,nn/rl^, (34) 

m n 

where A is the set of all charges to the left of the mid- 
plane, B is the set of all charges to the right of the mid- 
plane, Axjnn is the separation between the charges in 
the X direction, and rmn is the distance between the 
charges. In order to take into account lattice effects 
and correctly treat the periodic boundary conditions, we 
have calculated the force using the lattice coulomb force. 
This is done simply by replacing the continuum quan- 
tity Axmn/r'^n i^ Eq. |^ by the corresponding lattice 
expression 

4^y^ ^sin(2^fc^/L)e^"^^^"- 
4ELsin2(7rfc./L) 

where ki — 0, 1, 2, ....L — 1. 

Although the osmotic pressure could be calculated by 
observing the change in free energy as the separation of 
the plates is changed, work in Ref. [y has shown that the 
pressure calculated using Eq. has fewer fiuctuations 
than the pressure calculated using the free energy differ- 
ence. 

The pressure for univalent ions is shown in Fig. |2] for a 
range of plate separations. Results are shown from both 
the 50x50x50 lattice and the 100x50x50 lattice. The ions 
are moved using the coupled Metropolis update method 
described in Sec. IIIII as the heat bath method is diffi- 
cult to adapt to the asymmetric lattice. Figure 13 shows 
the first term of Eq. the concentration of ions on the 
plane between the plates, for the same simulations. The 
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Plate Separation (A) Plate Separation (A) 



FIG. 2: Osmotic pressure from simulations of univalent ions 
at a range of plate separations. 
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Plate Separation (A) 

FIG. 3: Concentration of charges at the plane between the 
charges from simulations of univalent ions. Left axis has the 
same units as the pressure results. Right axis gives the num- 
ber of charges on the center plane (50x50x50 lattice) or center 
two planes (100x50x50 lattice). This is one of the terms con- 
tributing to the osmotic pressure. 

left axis of this plot gives the average number of ions in 
a 1 A by 50 A by 50 A rectangular box centered between 
the plates. Figure0]shows the second term of Eq.|23l The 
differences in results from the two lattice sizes are mod- 
est, showing that the errors due to lattice discretization 
are small. Figure [S] shows the concentration profiles of 
the ions in solution from the simulations on the 50x50x50 
lattice. The ions are attracted to the plates, but a small 
density of ions remains in the center of the gap between 
the planes. 

For the divalent ions we only consider the 50x50x50 
lattice. Here we use the heat-bath method for mov- 
ing the particles and updating the electric fields. There 
are 10,000 equilibration steps and 200,000 measurement 



FIG. 4: Electrostatic force between the halves of the system 
for univalent ions. This is the second term contributing to 
the osmotic pressure. 
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FIG. 5: Concentration profile of the univalent ions in solution 
for a range of plate separations. For each x, the average 
number of ions in the system a distance x from the center 
plate is shown. The ion concentration is outside of the 
range shown. 



steps, and other parameters are the same as the univa- 
lent case. The solid line in Fig. shows the pressure in 
the divalent case, while dashed line and dashed-dotted 
line show the first and second terms of Eq. 123 Qual- 
itatively, our result agree with previous theoretical and 
Monte Carlo work, although a direct quantitative com- 
parison is difficult because of the differences in how the 
interaction with the plate is treated. The concentration 
profiles of the ions in solution for these simulations are 
shown in Fig.[7| Note that the ions are much more tightly 
bound to the plates in the divalent case as compared with 
the univalent case. 

The simulations discussed here would be difficult or 
impossible to perform with the simple particle Metropo- 
lis move, which does not adjust the field on neighboring 
plaquettes. By using the coupled Metropolis updating 
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FIG. 6: Results from simulations on divalent ions at a range 
of plate separations on a 50x50x50 lattice. Shown with a solid 
line is the total osmotic pressure divided by RT. Shown with 
a dashed line and a dotted-dashed line are the two terms that 
contribute to it, the ion concentration on the center plane and 
the electrostatic force between the two halves of the system. 
The axis on the left gives the values in units of micromolars; 
the axis on the right gives the concentration on the center 
plane in terms of the number of ions in our system on the 
center plane. 
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FIG. 7: Concentration profile of the divalent ions in solution 
for a range of plate separations. For each x, the average 
number of ions in the system a distance x from the center 
plate is shown. The ion concentration is outside of the 
range shown. 



move described in this work we are able to increase the 
acceptance rate for particle moves in the univalent case 
to 0.13 from 10~^, and in the divalent case to 0.0056 from 
less than 10~^. The anisotropy of the 100x50x50 lattice 
affects the acceptance rates. In this case the univalent 
acceptance rates using the coupled Metropolis move in- 
creases to 0.48 in the x direction, but the acceptance rate 
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FIG. 8: Comparison of pressure autocorrelation functions us- 
ing the coupled Metropolis move update, and the heat-bath 
particle move update. These autocorrelation functions are all 
from simulations with divalent ions. 



for particle moves perpendicular to the x direction drops 
to 0.0061. This asymmetry occurs because the coupled 
update procedure is no longer able to effectively spread 
the change in electric field when a particle is moved in 
the y or z directions. 

When using the heat-bath approach to moving the 
particles, the acceptance rate (0.10 for univalent ions 
and 0.0045 for divalent ions) is lower than the coupled 
Metropolis move acceptance rate, but still has a much 
greater acceptance rate than the simple Metropolis par- 
ticle move that does not adjust the electric field on neigh- 
boring plaquettes. The advantage of the heat-bath ap- 
proach is that it better decorrelates the system so that 
the observables have a shorter autocorrelation time. The 
autocorrelation function of a observable At is given by 



C{t) 



1 



N-t 



N ~t 



Y^{A,-A){A,+t-A), (36) 



where A is the average of At- The autocorrelation func- 
tion of the pressure is shown in Fig. |S| for plate sepa- 
rations of 8 and 14 with divalent ions using both the 
coupled Metropolis particle move update and the heat- 
bath particle move update. Although the autocorrela- 
tion time (obtained by integrating the autocorrelation 
function) increases with the larger plate separation, the 
autocorrelation time with the heat-bath update is con- 
sistently smaller than the autocorrelation time for the 
coupled Metropolis update. Across all plate separations, 
the autocorrelation times from heat-bath updates were 
five to ten times smaller than the autocorrelation times 
from coupled Metropolis updates. This more than com- 
pensates for the additional computational cost per sweep 
(approximately twice that of the coupled Metropolis up- 
date) of the heat-bath update. 
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V. CONCLUSIONS 

The development of efficient local algorithms for Monte 
Carlo simulation of Coulomb systems with non-uniform 
dielectric constants is crucial for the study of the larger 
and more physically realistic biophysical systems of in- 
terest to researchers 1^ ITollTT|. The technique of Maggs 
et al. shows great promise in fulfilling these goals, but 
it must be shown to be efficient and accurate in physi- 
cally interesting parameter ranges. Studying a system of 
parallel plates screened by ion with a large dimensionless 
inverse temperature, we see that the simplest method of 
moving particles, where the electric field is only modi- 



fied on the link traversed by the particle, gives unusably 
small acceptance rates. By updating the electric field 
on plaquettes neighboring the traversed link, we can in- 
crease the acceptance rates to a usable level. Using a 
heat bath approach reduces the autocorrelation time of 
the simulation. 
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